Holding water in a sieve—stable droplets without surface tension

Our understanding of supercritical fluids has seen exciting advances over the last decades, often in direct contradiction to established textbook knowledge. Rather than being structureless, we now know that distinct supercritical liquid and gaseous states can be distinguished and that a higher order phase transition - pseudo boiling - occurs between supercritical liquid and gaseous states across the Widom line. Observed droplets and sharp interfaces at supercritical pressures are interpreted as evidence of surface tension due to phase equilibria in mixtures, given the lack of a supercritical liquid-vapor phase equilibrium in pure fluids. However, here we introduce an alternative physical mechanism that unexpectedly causes a sharpening of interfacial density gradients in absence of surface tension: thermal gradient induced interfaces (TGIIF). We show from first principles and simulations that, unlike in gases or liquids, stable droplets, bubbles, and planar interfaces can exist without surface tension. These results challenge and generalize our understanding of what droplets and phase interfaces are, and uncover yet another unexpected behavior of supercritical fluids. TGIIF provide a new physical mechanism that could be used to tailor and optimize fuel injection or heat transfer processes in high-pressure power systems.

Imagine a droplet of ink suspended in water, and rather than diffusing outwards reducing the concentration gradients, the ink drop contracts and sharpens its interface instead! This is what we found can happen in the density field under certain conditions.
In this paper, we present thermodynamic conditions under which thermal gradient-induced interfaces (TGIIF) cause the density gradient to sharpen instead, without any surface force acting, seemingly in contradiction to our laws of diffusion. Furthermore, this interface sharpening occurs at supercritical pressures, in contradiction to van der Waals' gradient theory 1 and classical theories of nucleation 2 .
Injection at supercritical pressures is a highly relevant technology that is in widespread use in every Diesel engine 3 , every jet engine during take-off 4 , and in every main-stage rocket engine 5,6 . It is essential for new emerging high-efficiency carbon-capture power cycles [7][8][9] and new automotive combustion cycles 10,11 . Heat transfer to supercritical fluids is an essential technology in modern and projected power plants 12,13 , including solar 14 , nuclear 15 , and carbon capture 14 .
Over the last decade, it has become clear that a supercritical phase transition-pseudoboiling 16,17 -exists. Its main difference to subcritical boiling (or condensation) is the absence of an equilibrium coexistence of different phases; instead, the liquid-gas transition (or vice versa) occurs over a finite temperature interval 16,18,19 . Thus, despite early rejection of that idea 20 , even at supercritical pressures heat transfer deterioration can be caused by a (pseudo) boiling process [21][22][23] .
Droplet formation [24][25][26] at nominally supercritical pressures, on the other hand, has always been analyzed from a mixture thermodynamics standpoint. It is common for mixtures to exhibit a critical pressure that exceeds the pure fluid components' critical pressures 27 . Then, it is possible that a local mixture may find itself at a subcritical pressure, allowing phase separation, despite of the nominal supercritical pressure with respect to the pure components. This analysis has been performed for inert [28][29][30][31][32] and reactive mixing 33 . On a molecular scale, supercritical phase separation in mixtures can be attributed to the interfacial thickness 34 . All existing studies focus on mixtures [28][29][30][31][32]35,36 , where a phase equilibrium is assumed necessary for the formation of droplets.
However, the question arises as to whether the aforementioned pseudoboiling could provide a mechanism for droplets to form under true supercritical conditions, i.e., in absence of a phase equilibrium. In pseudoboiling, the liquid-gas transition occurs over a finite temperature interval across the pseudoboiling line 17 , an extension of the coexistence line. It is related to the Widom line [37][38][39][40] which is defined as locus of extrema in the thermodynamic response functions and thus can have many different characteristics depending on which response function is chosen 17 . The pseudoboiling line resolves this ambiguity with a precise definition based on the curvature of the free Gibbs enthalpy 17 . It is a fluid property and approximately marks the steepest isobaric thermal density gradient (∂ρ/∂T) p . In many injection problems, a spatial temperature gradient (∂T/∂x) exists between the droplet and its surroundings [24][25][26][28][29][30][31][32] . Thus, the pseudoboiling line seems a likely candidate to induce a maximum in the spatial density gradient (∂ρ/∂x), i.e., a density inflection point, in the presence of a spatial temperature gradient-much like what we see in subcritical droplets 34 .
This paper introduces and analyzes an alternative mechanism that stabilizes droplets and bubbles at supercritical pressures in the absence of surface tension, based on heat transfer rather than surface forces: thermal gradient-induced interfaces (TGIIF) Specifically, we will show that, even in a pure fluid, a temperature difference between a cold/ dense region and its warmer/lighter environment can be enough to cause a self-steepening and self-stabilizing density gradient that would be indistinguishable from acting surface tension in a shadowgraph.

Results
A simplified view of a phase interface is a discontinuous transition between two fluid properties. Figure 1a illustrates this view, where at some location ℓ an interface marks the switch from a high density to a low density. More physically and accurately, however, is the gradual transition seen in a real vapor-liquid interface 41,42 where the position of the interface can be found in the vicinity of the density inflection point, in contrast to the density profile of a gaseous diffusion process.
Definition 1 (Interface). Thus, here we consider an interface the location at which the density distribution exhibits a spatial inflection point.
Definition 2 (Stable). We consider stability the tendency of some system to approach or revert to some preferred state when it is perturbed from it.
We will now demonstrate that both properties can be simultaneously present in fluids without surface tension. a While an idealized liquid-vapor interface is often considered discontinuous, the real interface is smooth and characterized by a density inflection point 41 . b With no interface, e.g., in a gas, density in the presence of a temperature gradient will drop without an inflection point. c Equation (9) identifies states in which a density inflection point can be caused by a temperature gradient. CL is the coexistence line, T pb is the pseudoboiling temperature, T * is the temperature at which the thermal conductivity reaches a minimum, T α is the temperature of maximum thermal expansion. Middle row: Examples of fluid properties for d gas, e transcritical, and f very high-pressure fluid. The shaded area fulfils Eq. (9). Bottom: corresponding density fields at the respective conditions. Unlike the gas in (g), the transcritical condition (h) exhibits a very clear density interface. Surprisingly, even at very high pressures (i) where the fluid density no longer exhibits an inflection point with increasing temperature, a broad but sustained interface establishes.

Existence of a steady interface
We start with the analysis of the steady interface at constant pressure. When the interface is characterized by an inflection point in the spatial density distribution, it is linked to the thermodynamic thermal density gradient via the spatial temperature gradient, Then, using the chain rule for the 2nd derivative, setting it to zero, and using the convention of subscripts denoting differentiation, such that ρ T = (∂ρ/∂T) and ρ TT = (∂ 2 ρ/∂T 2 ), yields Using Fourier's law of heat conduction in 1D to eliminate T x results in the steady-state relation Equation (5) needs to be fulfilled for a density inflection point, and thus interface, to exist. The RHS describes heat transfer characteristics, whereas the LHS represents fluid thermodynamic properties. The simplest way to analyze Eq. (5) is to check for sign consistency. The squared term on the RHS is always positive. This means that an inflection point for fluids that expand when heated, i.e., ρ T < 0, can exist if and only if Further, through differentiation of Eq. (4) for steady state with constant q ″ throughout the field, where k > 0. Without loss of generality, we regard the case of Fig. 1a with a positive temperature gradient T x > 0 and find Finally, with k x = (∂k/∂T)(∂T/∂x), we obtain as a necessary condition for the existence of an interface. Note that Eq. (9) is a criterion based entirely on fluid properties, and does not depend on experimental boundary conditions. Figure 1c shows the region in which Eq. (9) is fulfilled in a pressure-temperature diagram.

Fluid state case evaluations
Using Eq. (9), we can now evaluate different fluid states for their potential to develop interfaces, results are shown in Fig. 1. For simplicity, we will restrict ourselves to isobaric cases, it is instructive to view the fluid property diagrams in "Methods". Here, T pb is the pseudoboiling temperature 16 which characterizes the supercritical liquid-gas transition and is defined as the locus of the isobaric heat capacity peaks. k min is the minimum value of the thermal conductivity evaluated at constant pressure and T * is the temperature at which k min is reached.
In an ideal gas with ρ = p=RT, sgnðρ TT Þ > 0, kinetic theory suggests a dependence k / ffiffiffi ffi T p , thus sgnðk T Þ > 0. As this contradicts Eq. (9), an ideal gas cannot exhibit an interface.
The same analysis can be used to show that interfaces cannot occur in liquids with typically sgnðρ TT Þ < 0 and sgnðk T Þ < 0.
In a transcritical fluid, for T < T pb , sgnðρ TT Þ < 0 and sgnðk T Þ < 0 like in a liquid, preventing an interface. On the other hand, for T > T * , with T * :¼ Tðk min Þ = T * , sgnðρ TT Þ < 0 and sgnðk T Þ > 0, like in a gas, likewise preventing an interface. This leaves T pb < T < T * as a region of unique fluid properties with sgnðρ TT Þ > 0 and sgnðk T Þ < 0, allowing for an interface to form. Figure 1f shows in a p-T diagram that an interface can exist at high subcritical and even very high supercritical pressures. In fact, we did not find an upper pressure limit to the existence of a density inflection point.

Analytical analysis and steady 1D simulation
A supercritical steady interface can be analyzed analytically. Consider two parallel walls at different temperatures, a cold T c and a hot T h . The space between them is filled with a pure fluid and the pressure is kept constant. The fluid does not flow, so heat is only transferred by conduction. This will allow our heat transfer to be analyzed by Fourier's law, Eq. (4). In fact, the temperature distribution is solely governed by the boundary temperatures and the fluid thermal conductivity.
An accurate modeling of transcritical fluid properties is an ongoing problem 43,44 . Here, we approximate the complex thermal conductivity distribution shown in Fig. 1 as piece-wise linear, with a liquid branch and a gaseous branch, a detailed derivation can be found in the methodology. Then, exact analytical solutions to the temperature distribution can be determined from Fourier's law, Somewhat unexpectedly, agreement between the analytical and the numerical solutions is excellent for a wide range of conditions. Beginning with the ideal gas, the linear model closely matches the density distribution of the simulation, which exhibits no inflection point and thus no interface. At transcritical conditions (7 MPa, p/p cr ≈ 1.2), a very distinct inflection point develops, separating the domain into a liquid-filled and a gas-filled side. Even at extremely high pressures of 100 MPa (p/p cr ≈ 20) does the inflection point not vanish; however, the density gradient has become very smooth and a clear separation into liquid and gaseous has all but vanished.

Transient gradient sharpening
To study the transient interface, we perform a 2D parametric study looking at droplets of n-heptane at several different pressures and environmental temperatures.  Figure 2 demonstrates that TGIIF are indeed acting to create stable interfaces: Fig. 2a illustrates how an initially fuzzy droplet develops a sharper interface over time; Fig. 2b demonstrates the regression of an elliptical cross-section toward a circular cross-section; finally, Fig. 2c shows a supercritical droplet without TGIIF, diffusing, as one would expect for all of them.
In order to analyze the vaporization transient, we determine a representative timescale from dimensional analysis of transient heat conduction, c.f. "Time scale analysis", The associated timescale is then We define the nondimensional time τ where α is the minimum thermal diffusivity for a fluid at a given pressure, θ is the temperature at which the minimum thermal diffusivity occurs, ΔT is the temperature between the warm gas and the cold liquid, and D is the initial diameter of the droplet. The transient development of the maximum density gradient is quantified in Fig. 2a and b for two different pressures. We can clearly see how the density gradient grows linearly until a maximum is reached, if the temperature difference ΔT between droplet and surroundings is large enough. Figure 2f maps the positions at which TGIIF and diffuse interfaces are found.

Transient vaporization
The dimensional analysis in Eq. (13) suggests that supercritical vaporization will be linear in D 2 , akin to Spalding's D 2 law in subcritical droplet evaporation 46 . Following Fig. 2, we consider the interface of the droplet the position of the maximum density gradient, then D is the diameter based on that location. Figure 3a-c shows how the square of the diameter D, nondimensionalized with the initial diameter D 0 , indeed shrinks linearly in time for TGIIF conditions-much like a subcritical droplet! For sufficiently high ΔT, the curves collapse; for T r < 2 and p r < 4, the droplet grows initially before following the D 2 relation. Figure 3d, e reveals another unexpected similarity of TGIIF to subcritical droplet evaporation: Regardless of the environmental temperature, the temperature at the interface matches a single temperature-the pseudoboiling temperature T pb , which can be considered the supercritical complement to subcritical saturation temperature, and is found on the pseudoboiling line as an extension to the coexistence line 16,17 . Pseudoboiling seizes to be a relevant phenomenon for p r > 3, and Fig. 3f shows a larger deviation from T pb , although it still appears to act as a (weak) attractor.
We now turn to the droplet lifetimes as our final analysis. Figure. 3g, h illustrates the droplet evolution using snapshots at different stages of vaporization, for a stable TGIIF droplet and a diffuse droplet, respectively. We consider the droplet lifetime the time until the maximum density gradient has reached its minimum value or a zero radius. We can then compile and compare droplet lifetimes over the whole parameter range, i.e., T r = [1, 1.5, 2, 2.5], p r = [1.5, 2, 4, 6], D = [0.001, 0.01, 0.1] for all cases that exhibited TGIIF.
A dimensional analysis, such as the one performed to obtain Eq. (13) and as presented in more detail in the methods part, can reveal relations only up to a proportionality constant. We thus introduce a vaporization number Va as the nondimensional proportionality constant to Eq. (13), We find that for Va = 0.31 the predicted droplet lifetime is accurate to within 50% over five decades of lifetimes obtained from simulations. Figure 3i illustrates this. Thus, Va = 0.31 seems to be the characteristic vaporization number for the problem at hand.

Discussion
We introduce a new physical mechanism that acts to stabilize liquid-gas interfaces, such as droplets or bubbles, in the absence of phase equilibrium and surface tension: thermal gradient-induced interfaces (TGIIF). TGIIF are surprisingly similar to classical droplet evaporation: TGIIF droplets exhibit sharpening density gradients and cause elliptical cross-sections to regress towards circular cross-sections; they vaporize following a D 2 law, while their interfaces hold the supercritical pseudoboiling temperature, akin to the subcritical saturation temperature.
Our first principle analysis of criteria for a stable density inflection point yields as a necessary condition that the concavity of the thermal density gradient and the slope of the thermal conductivity have opposite signs, sgnðρ TT Þ = À sgnðk T Þ. This condition holds under transcritical conditions; it does not hold in liquids or gases.
The nondimensional vaporization number Va relates droplet size, temperature gradients, and thermal diffusivity to droplet lifetime. We find that Va = 0.31 matches simulated droplet lifetimes over five decades, within 50% error.
Ultimately, this work questions and extends the definition of what constitutes a droplet or a bubble. Photographic evidence of stable spherical structures in high-pressure injection does not necessarily imply the existence of a phase equilibrium or surface tension. Our study answers the long standing question as to whether transcritical interfaces can be stable: they can, even in absence of surface tension. It provides justification for two-layer models for supercritical heat transfer analysis.

Solver
The CFD solver used was open-source SU2 47,48 , more specifically the low-mach approximation solver was used. The convection scheme used is the flux difference splitting scheme (FDS). The FDS scheme is an upwind scheme and typically has first-order accuracy, but second-order accuracy is achieved via Monotonic Upwind Scheme for Conservation Laws (MUSCL) 47 . The solver was then extended with tiny neural networks (TNN) in order to model efficiently and accurately capture the complex supercritical fluid properties. An in-depth explanation of the TNN and its implementation into SU2 are found in refs. 22,45, representative fluid property distributions are shown for n-heptane and oxygen in Fig. 4.
The case was set up so the left half of the domain was the cold liquid side with cold wall of the same temperature on the edge, and on the right side was the warm vapor with a hot wall of the same temperature on the edge. The heat would then diffuse from the warm side to the cold side until a constant heat flux was reached across the domain and that was the steady state.

Analytical model methodology
A supercritical droplet or bubble structure should be described analytically. To portray a supercritical droplet, Fourier's Law is used across a simple one-dimensional steady-state case with a linear approximation of the fluid thermal conductivity. The linear model allows for hand derivations of an equation to describe the temperature profile which is then coded in Python to iterate across a spatial domain and plot the steady-state temperature profile.
Deriving a pseudoboiling distance. The analytical case is defined as a cold wall of temperature T c and a hot wall of temperature T h separated by a distance L. The steady-state heat flux is governed by Fourier's Law.
For a 1D model, the gradient is equal to d T/d x. Thermal conductivity is not constant and can be expressed as a function of temperature.
Àq 00 = kðTÞ dT dx ð16Þ Next, a fluid will be chosen at a pressure which pseudoboiling occurs. Temperatures T c and T h are then defined as below and above the pseudoboiling temperature T pb , respectively. This is intended to capture the pseudoboiling transition. A linear function is then fit to the thermal conductivity such that k c (T) spans from T c to some arbitrary interface temperature,T. This region will be referred to as the cold side. Also, a linear function k h (T) approximates the region fromT to T h , the hot side. To create a linear equation for each region, a temperature between the wall and interface temperatures is chosen T region for both the hot and cold sides. These line-fitting temperatures have a corresponding thermal conductivity, k region .
and b region = k wall À m region T wall ð19Þ In this model, there exists anx whereT occurs. Because the model is in a steady state, the heat flux is constant across any region d x. So, the heat flux across the region below the transition temperature shall equal the heat flux in the region above the transition temperature and k(T) can be substituted in Next, both sides are integrated over respective temperatures and distances, and the integral is solved The equations can be rearranged and the large ratio of integrated k(T) d T represented as ξ.
Rearranged with ξ And rearranged again to isolatex With (25),x is known for a fully defined case with: constant heat flux; a particular isobaric fluid; hot, cold, and pseudoboiling temperatures; and a discrete region x c to x h .
Plotting a temperature profile. To create a temperature profile, Fourier's Law (16) can be integrated across the domain, 0 to x, and from T c to T(x) (26). Because k(T) is different across the two regions, the results will differ depending on what value of x is chosen.
For x <=x the integration results in (27) where T(x) can be solved for with the quadratic formula, and positive results chosen in (28).
For x >x, we can integrate over the hot region, which yields in a similar formula With an expression for T(x), the spatial temperature gradient can be found as well by finding the derivative of T(x) in respect to x. This yields the following two equations.
1D steady state The analytical model outlined in "Analytical model methodology" can be applied to oxygen at 7 MPa between 100 K and 300 K walls, across a unit domain of 1 m to yield fractional results. For comparison to a computational solver, SU2, the relative domain results are applied to a 1-mm distance between hot and cold walls.
To best represent the thermal conductivity of 7-MPa oxygen, the equations for the thermal conductivity model (k-model) were chosen to reflect the liquid region between 100 K and 158 K for the cold side, and between 250 K and 300 K to reflect the gaseous region. The resulting k-model is shown below in Fig. 5 and a sample of cases with parameters of interest can be found in Table 1.

2D droplets and bubbles
Mesh. In Figure 6 are pictures of the mesh used for the 2D droplet simulations. The left picture shows the whole domain and the right picture zooms in on the center of the mesh to show the refined part of the mesh where the interface occurs.

Time-scale analysis
We identify a representative timescale of the vaporization process by performing dimensional analysis. 1D heat conduction through a constant property medium is governed by We define reference values and corresponding nondimensional variables such that x = x=D and t = t=τ. For the temperature term on the lhs we introduce a reference temperature T = T=θ, on the rhs we choose a reference temperature difference instead, T = T=ΔT. Then, Article https://doi.org/10.1038/s41467-023-39211-z Eq. (32) can be rewritten and we identify a nondimensional group which we will dub the vaporization number The associated timescale is then τ is used in our analysis to nondimensionalize the time. α is the minimum thermal diffusivity for a fluid at a given pressure, θ is the temperature at which the minimum thermal diffusivity occurs, ΔT is the temperature between the warm gas and the cold liquid, and D is the initial diameter of the droplet.

Data availability
The fluid property data used for this study was through Coolprop 49 and NIST 50 , and the data is available through their websites. The data can be requested from the authors.

Code availability
The tiny neural networks used for fluid property modeling is the only custom code that was used and it is available through the GitLab website https://gitlab.com/tfxlab/thermoml.git. The solver used for  the study, open-source SU2, is available on their GitHub website https://github.com/su2code/SU2. Article https://doi.org/10.1038/s41467-023-39211-z